{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Traveling Salesman Problem\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "In this example, you’ll learn how to tackle one of the most famous combinatorial optimization problems in existence: the Traveling Salesman Problem (TSP). The goal of the TSP – to find the shortest possible route that visits each city once and returns to the original city – is simple, but solving the problem is a complex and challenging endeavor. We’ll show you how to do it!\n",
    "\n",
    "This modeling example is at the advanced level, where we assume that you know Python and the Gurobi Python API and that you have advanced knowledge of building mathematical optimization models. Typically, the objective function and/or constraints of these examples are complex or require advanced features of the Gurobi Python API.\n",
    "\n",
    "**Download the Repository** <br />\n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Motivation\n",
    "\n",
    "The Traveling Salesman Problem (TSP) is one of the most famous combinatorial optimization problems. This problem is very easy to explain, but very complicated to solve – even for instances with a small number of cities. More detailed information on the TSP can be found in the book The Traveling Salesman Problem: A Computational Study [1], or at the TSP home page [2]. If you are interested in the history and mathematical background of the TSP, we recommend that you watch the video by William Cook [3].\n",
    "\n",
    "The origin of the traveling salesman problem is not very clear; it is mentioned in an 1832 manual for traveling salesman, which included example tours of 45 German cities but was not formulated as a mathematical problem. However, in the 1800s, mathematicians William Rowan Hamilton and Thomas Kirkman devised mathematical formulations of the problem.\n",
    "\n",
    "It seems that the general form of the Traveling Salesman Problem was first studied by Karl Menger in Vienna and Harvard in the 1930s.\n",
    "\n",
    "The problem became more and more popular in the 1950s and 1960s. In particular, George Dantzig, D. Ray Fulkerson, and Selmer M. Johnson at the RAND Corporation solved the 48-state problem by formulating it as a linear programming problem. The methods they described in their paper on this topic set the foundation for future work in combinatorial optimization, especially highlighting the importance of cutting planes.\n",
    "\n",
    "In the early 1970s, the concept of P vs. NP problems created excitement in the theoretical computer science community. In 1972, Richard Karp demonstrated that the Hamiltonian cycle problem was NP-complete, implying that the traveling salesman problem was NP-hard.\n",
    "\n",
    "Increasingly sophisticated codes led to rapid increases in the sizes of the traveling salesman problems solved. Dantzig, Fulkerson, and Johnson had solved a 48-city instance of the problem in 1954. Martin Grötschel more than doubled this 23 years later, solving a 120-city instance in 1977. Harlan Crowder and Manfred W. Padberg again more than doubled this in just 3 years, with a 318-city solution.\n",
    "\n",
    "In 1987, rapid improvements were made, culminating in a 2,392-city solution by Padberg and Giovanni Rinaldi. In the following two decades, great strides were made with David L. Applegate, Robert E. Bixby, Vasek Chvátal, & William J. Cook solving a 3,308-city instance in 1992, a 7,397-city instance in 1994, a 24,978-city instance in 2004, and an 85,900-city instance in 2006 – which is the largest 2-D Euclidean TSP instance ever solved. William Cook et. al. wrote a program called Concorde TSP Solver for solving the TSP [4]. Concorde is a computer code for the symmetric TSP and some related network optimization problems. The code is written in the ANSI C programming language and it has been used to obtain the optimal solutions to the full set of 110 TSPLIB instances, the largest instance is a 109,399 node 3-D “star” instance.\n",
    "\n",
    "The continued interest in the TSP can be explained by its success as a general engine of discovery and a steady stream of new applications. Some of the general applications of the TSP are as follows:\n",
    "* Scheduling and routing problems.\n",
    "* Genome sequencing.\n",
    "* Drilling problems.\n",
    "* Aiming telescopes and x-rays.\n",
    "* Data clustering.\n",
    "* Machine scheduling.\n",
    "\n",
    "We use this classic combinatorial optimization problem to demonstrate how Gurobi can be used to easily and effectively solve small-sized problem instances of the TSP. However, in order to be able to solve larger instances, one needs more sophisticated techniques – such as those implemented in the Concord TSP Solver."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Description\n",
    "The TSP can be defined as follows: for a given list of cities and the distances between each pair of them, we want to find the shortest possible route that goes to each city once and returns to the origin city.\n",
    "\n",
    "There is a class of Traveling Salesman Problems that assumes that the distance of going from city $i$ to city $j$  is the same as going form city $j$ to city $i$, this type of Traveling Salesman Problem  is also known as the symmetric Traveling Salesman Problem. In this example, we use Euclidean distances, but the TSP model formulation is valid independent of the way in which the individual distances are determined.\n",
    "\n",
    "\n",
    "## Solution Approach\n",
    "\n",
    "Mathematical programming is a declarative approach where the modeler formulates a mathematical optimization model that captures the key aspects of a complex decision problem. The Gurobi Optimizer solves such models using state-of-the-art mathematics and computer science.\n",
    "\n",
    "A mathematical optimization model has five components, namely:\n",
    "\n",
    "* Sets and indices.\n",
    "* Parameters.\n",
    "* Decision variables.\n",
    "* Objective function(s).\n",
    "* Constraints.\n",
    "\n",
    "We now present a MIP formulation of the TSP that identifies the shortest route that goes to all the cities once and returns to the origin city.\n",
    "\n",
    "## TSP Model Formulation\n",
    "\n",
    "### Sets and Indices\n",
    "$i, j \\in Capitals $: indices and set of US capital cities.\n",
    "\n",
    "$\\text{Pairings}= \\{(i,j) \\in Capitals \\times Capitals \\}$: Set of allowed pairings\n",
    "\n",
    "$S \\subset Capitals$: A subset of the set of US capital cities.\n",
    "\n",
    "$G = (Capitals, Pairings)$: A graph where the set $Capitals$ defines the set of nodes and the set $Pairings$ defines the set of edges. \n",
    "\n",
    "### Parameters \n",
    "\n",
    "$d_{i, j} \\in \\mathbb{R}^+$: Distance from capital city $i$ to capital city $j$, for all $(i, j) \\in Pairings$. \n",
    "\n",
    "Notice that the distance from capital city $i$ to capital city $j$ is the same as the distance from capital city $j$ to capital city $i$, i.e. $d_{i, j} = d_{j, i}$. For this reason, this TSP is also called the symmetric Traveling Salesman Problem.\n",
    "\n",
    "### Decision Variables\n",
    "$x_{i, j} \\in \\{0, 1\\}$: This variable is equal to 1, if we decide to connect city $i$ with city $j$. Otherwise, the decision variable is equal to zero.\n",
    "\n",
    "### Objective Function\n",
    "- **Shortest Route**. Minimize the total distance of a route. A route is a sequence of capital cities where the salesperson visits each city only once and returns to the starting capital city.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Min} \\quad Z = \\sum_{(i,j) \\in \\text{Pairings}}d_{i,j} \\cdot x_{i,j}\n",
    "\\tag{0}\n",
    "\\end{equation}\n",
    "\n",
    "### Constraints \n",
    "- **Symmetry Constraints**. For each edge $(i,j)$, ensure that the city capitals $i$ and $j$ are connected, if the former is visited immediately before or after visiting the latter.\n",
    "\n",
    "\\begin{equation}\n",
    "x_{i, j} = x_{j, i} \\quad \\forall (i, j) \\in Pairings\n",
    "\\tag{1}\n",
    "\\end{equation}\n",
    "\n",
    "- **Entering and leaving a capital city**. For each capital city $i$, ensure that this city is connected to two other cities. \n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{(i,j) \\in \\text{Pairings}}x_{i,j} = 2 \\quad \\forall  i \\in Capitals\n",
    "\\tag{2}\n",
    "\\end{equation}\n",
    "\n",
    "- **Subtour elimination**. These constraints ensure that for any subset of cities $S$ of the set of $Capitals$, there is no cycle. That is, there is no route that visits all the cities in the subset and returns to the origin city.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{(i \\neq j) \\in S}x_{i,j} \\leq |S|-1 \\quad \\forall  S \\subset  Capitals\n",
    "\\tag{3}\n",
    "\\end{equation}\n",
    "\n",
    "- **Remark**. In general, if the number of cities of the TSP is $n$, then the possible number of routes is n\\!.\n",
    "Since there are an exponential number of constraints ($2^{n} - 2$) to eliminate cycles, we use lazy constraints to dynamically eliminate those cycles. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Implementation\n",
    "\n",
    "Consider a salesperson that needs to visit customers at each state capital of the continental US. The salesperson wants to identify the shortest route that goes to all the state capitals.\n",
    "\n",
    "This modeling example requires the following libraries that are not part of the standard Python distribution:\n",
    "* **folium**: to create maps.\n",
    "* **gurobipy**: provides Gurobi algorithms to solve MIP models.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading Input Data\n",
    "The capital names and coordinates are read from a json file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# Read capital names and coordinates from json file\n",
    "try:\n",
    "  capitals_json = json.load(open('capitals.json'))\n",
    "# when running locally the following lines can be omitted\n",
    "except: \n",
    "  import urllib.request\n",
    "  url = 'https://raw.githubusercontent.com/Gurobi/modeling-examples/master/traveling_salesman/capitals.json'\n",
    "  data = urllib.request.urlopen(url).read()\n",
    "  capitals_json = json.loads(data)\n",
    "\n",
    "capitals = []\n",
    "coordinates = {}\n",
    "for state in capitals_json:\n",
    "    if state not in ['AK', 'HI']:\n",
    "      capital = capitals_json[state]['capital']\n",
    "      capitals.append(capital)\n",
    "      coordinates[capital] = (float(capitals_json[state]['lat']), float(capitals_json[state]['long']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data computation\n",
    "The following function calculates the distance for each pair of state capitals. Since we are solving the _symmetric_ traveling salesman problem, we use _combinations_ of cities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from itertools import combinations\n",
    "\n",
    "# Compute pairwise distance matrix\n",
    "\n",
    "def distance(city1, city2):\n",
    "    c1 = coordinates[city1]\n",
    "    c2 = coordinates[city2]\n",
    "    diff = (c1[0]-c2[0], c1[1]-c2[1])\n",
    "    return math.sqrt(diff[0]*diff[0]+diff[1]*diff[1])\n",
    "\n",
    "dist = {(c1, c2): distance(c1, c2) for c1, c2 in combinations(capitals, 2)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Code\n",
    "We now write the model for the TSP, by defining decision variables, constraints, and objective function. Because this is the _symmetric_ traveling salesman problem, we can make it more efficient by setting the _object_ x[j,i] to x[i,j], instead of a constraint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Python 3.7 & Gurobi 9.0.0\n",
    "\n",
    "m = gp.Model()\n",
    "\n",
    "# Variables: is city 'i' adjacent to city 'j' on the tour?\n",
    "vars = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='x')\n",
    "\n",
    "# Symmetric direction: use dict.update to alias variable with new key\n",
    "vars.update({(j,i):vars[i,j] for i,j in vars.keys()})\n",
    "\n",
    "# Constraints: two edges incident to each city\n",
    "cons = m.addConstrs(vars.sum(c, '*') == 2 for c in capitals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Callback Definition\n",
    "Subtour constraints prevent multiple loops in a TSP tour. Because there are an exponential number of these constraints, we don't want to add them all to the model. Instead, we use a callback function to find violated subtour constraints and add them to the model as lazy constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Callback - use lazy constraints to eliminate sub-tours\n",
    "\n",
    "def subtourelim(model, where):\n",
    "    if where == GRB.Callback.MIPSOL:\n",
    "        # make a list of edges selected in the solution\n",
    "        vals = model.cbGetSolution(model._vars)\n",
    "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
    "                             if vals[i, j] > 0.5)\n",
    "        # find the shortest cycle in the selected edge list\n",
    "        tour = subtour(selected)\n",
    "        if len(tour) < len(capitals):\n",
    "            # add subtour elimination constr. for every pair of cities in subtour\n",
    "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
    "                         <= len(tour)-1)\n",
    "\n",
    "# Given a tuplelist of edges, find the shortest subtour\n",
    "\n",
    "def subtour(edges):\n",
    "    unvisited = capitals[:]\n",
    "    cycle = capitals[:] # Dummy - guaranteed to be replaced\n",
    "    while unvisited:  # true if list is non-empty\n",
    "        thiscycle = []\n",
    "        neighbors = unvisited\n",
    "        while neighbors:\n",
    "            current = neighbors[0]\n",
    "            thiscycle.append(current)\n",
    "            unvisited.remove(current)\n",
    "            neighbors = [j for i, j in edges.select(current, '*')\n",
    "                         if j in unvisited]\n",
    "        if len(thiscycle) <= len(cycle):\n",
    "            cycle = thiscycle # New shortest subtour\n",
    "    return cycle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m._vars = vars\n",
    "m.Params.lazyConstraints = 1\n",
    "m.optimize(subtourelim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "We retrieve the optimal solution of the TSP and verify that the optimal route (or tour) goes to all the cities and returns to the origin city."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve solution\n",
    "\n",
    "vals = m.getAttr('x', vars)\n",
    "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
    "\n",
    "tour = subtour(selected)\n",
    "assert len(tour) == len(capitals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimal route is displayed in the following map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Map the solution\n",
    "\n",
    "import folium\n",
    "\n",
    "map = folium.Map(location=[40,-95], zoom_start = 4)\n",
    "\n",
    "points = []\n",
    "for city in tour:\n",
    "  points.append(coordinates[city])\n",
    "points.append(points[0])\n",
    "\n",
    "folium.PolyLine(points).add_to(map)\n",
    "\n",
    "map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.dispose()\n",
    "gp.disposeDefaultEnv()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "The Traveling Salesman Problem (TSP) is the most popular combinatorial optimization problem. This problem is very easy to explain, although it is very complicated to solve. The largest TSP problem solved has 85,900 cities. The TSP is a source of discovery for new approaches to solve complex combinatorial optimization problems and has led to many applications.\n",
    "\n",
    "In this modeling example, we have shown how to formulate the symmetric Traveling Salesman Problem as a MIP problem. We also showed how to dynamically eliminate subtours by using lazy constraints. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] D. L. Applegate, R. E. Bixby, V. Chvatal and W. J. Cook , The Traveling Salesman Problem: A Computational Study, Princeton University Press, Princeton, 2006.\n",
    "\n",
    "[2] http://www.math.uwaterloo.ca/tsp/index.html\n",
    "\n",
    "[3] https://www.youtube.com/watch?v=q8nQTNvCrjE&t=35s\n",
    "\n",
    "[4] http://www.math.uwaterloo.ca/tsp/concorde.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
